R Notebook

Initiating

library(oro.dicom)
library(tidyverse)
library(plotly)
library(data.table)
library(htmlwidgets)
library(coreCT)
library(magick)
library(raster)
library(doParallel)
library(oro.dicom)
library(magick)
library(imager)
ref<-1000*data.frame(hmax=(4.64 +0.77),hmin=(4.64 -0.77), mmax=(3.14 + .65), mmin =(3.14- .65 ))

osis<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/osic-pulmonary-fibrosis-progression")
osistrain<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/osic-pulmonary-fibrosis-progression/train")
osistest<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/osic-pulmonary-fibrosis-progression/test")
osisimg<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/img")

setwd(osis)

#invalids<-"ID00011637202177653955184"   #https://www.kaggle.com/rashmibanthia/corrupt-files

trd<-read.csv("train.csv")  %>% dplyr::rename(Week=Weeks)
#subset(Patient!= invalids)
ted<-read.csv("test.csv") %>% dplyr::rename(Week=Weeks)
invalids<-c("ID00011637202177653955184","ID00052637202186188008618", "ID00128637202219474716089","ID00132637202222178761324")

dd<-list.files(osistrain)
dd <- setdiff(dd,invalids)

*** load dicom img and convert to HU units

x=dd[1]

There were 32 warnings (use warnings() to see them)

read_x<-function(x,thr){
  path<-file.path(osistrain,x)
    s<-readDICOM(path)
    #img<-s$img[[index]] 
    
    ct.slope <- unique(extractHeader(s$hdr, "RescaleSlope"))
    ct.int   <- unique(extractHeader(s$hdr, "RescaleIntercept")) 
    Ct.n<- extractHeader(s$hdr, "InstanceNumber")
    ss <- lapply(s$img, function(x) x*ct.slope + ct.int) 
   # ff<-function(x){x [ x < -2000 ]<- NA;return(x)}
    # trim background
    ss<-  lapply(ss,function(x){x[x< thr]<-F;return(x)})
    #reorder:
    fo<-function(i){
      index<-which(Ct.n ==i)
      return(ss[[index]])
    }
    so<-lapply(1:length(ss),fo)
    return(so)
}
so<-read_x(x,-2000)
    

**imager

l.im<-lapply(t(so[1:3]),as.cimg)
img<-as.imlist(l.im) ; names(img)<-NULL
plot(img,axes = F)

img<-t(so[[7]])



im<-as.cimg( img )
plot(im)

*** more imager processing

px<-(im < - 600 & im > -860 )

plot(px)

plot(im)
highlight(px)

NA
NA

plot(im)
px.flood(im,78,296,sigma=.21) %>% highlight

pal <- colorRamp(c("black","white"))

fig <- plot_ly(z = ~ so[[7]], colors = pal) %>%
add_surface(
  contours = list(
    z = list(
      show=TRUE,
      usecolormap=TRUE,
      highlightcolor="#ff0000",
      project=list(z=TRUE)
    )
  )
)
fig <- fig %>% layout(
  scene = list(
    camera=list(
      eye = list(x=1.87, y=0.88, z=-0.64)
    )
  )
)

fig

NA

** new threshould to focus on lungs

sfail<-read_x(x,-600)

'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2

pal <- colorRamp(c("black","white"))

fig <- plot_ly(z = ~ sfail[[7]], colors = pal) %>%
add_surface(
  contours = list(
    z = list(
      show=TRUE,
      usecolormap=TRUE,
      highlightcolor="#ff0000",
      project=list(z=TRUE)
    )
  )
)
fig <- fig %>% layout(
  scene = list(
    camera=list(
      eye = list(x=1.87, y=0.88, z=-0.64)
    )
  )
)

fig

It has cut most of the lungs…that is below -800

**** lets start lung segmentation linear model:

d <- as.data.frame(im)
##Subsamble, fit a linear model
m <- sample_n(d,1e4) %>% lm(value ~ x*y,data=.) 
##Correct by removing the trend
im.c <- im-predict(m,d)
out <- threshold(im.c)
plot(out)

 out<-threshold(im,"auto")
out <- clean(out,3) %>% imager::fill(7)
plot(im)
highlight(out)

NA
NA

im[ out]<- NA
plot(im)

 plot_ly(
            z = so[[7]],
            #colorscale = list(c(0,0.5,1),c("blue", "white", "red")),
            colors=pal,
            type = "heatmapgl"
            )

NA
NA

a<-kmeans(out,2)
colors<- a$centers[a$cluster[1,],]

im<-as.cimg(so[[7]]) %>% imrotate(-90)


detect.edges <- function(im,sigma=1)
    {
        isoblur(im,sigma) %>% imgradient("xy") %>% enorm %>% imsplit("c") %>% add
    }

edges <- detect.edges(im,2) %>% sqrt 
plot(edges)

pmap <- 1/(1+edges) #Priority inv. proportional to gradient magnitude
plot(pmap,main="Priority map") #Nice metal plate effect! 

NA
NA

seeds <- imfill(dim=dim(pmap)) #Empty image
seeds[170,233,1,1] <- 1 #Background pixel 
seeds[36,255,1,1] <- 2 #Foreground pixel

wt1 <- watershed(seeds,pmap)
plot(wt1,main="Watershed segmentation")

the other lung

seeds <- imfill(dim=dim(pmap)) #Empty image
seeds[355,247,1,1] <- 1 #Background pixel 
seeds[56,252,1,1] <- 2 #Foreground pixel

wt2 <- watershed(seeds,pmap)
plot(wt2,main="Watershed segmentation")

NA
NA

now, both

wt <- wt1 + wt2

plot(wt)

NA
NA

wt[wt==3]<-T
wt[wt==4]<-F
lungs<-wt * im
plot(lungs)

NA
NA

dx <- imgradient(im,"x")
dy <- imgradient(im,"y")
grad.mag <- sqrt(dx^2+dy^2)
plot(grad.mag,main="Gradient magnitude")

** lets try kmeans

a<-kmeans(im)
plot(im) %>% highlight(a$centers)